## ============================================================================
## 10_si_s2_practitioner_heterogeneity.R
##
## Purpose:   Treatment effect heterogeneity analysis for Study 2 using
##            Generalized Random Forests (causal forests).
##            Corresponds to SI Section S2.3.1.
##
## Inputs:    combined_data.rds (via 00_setup.R)
## Outputs:   grf_calibration_au_study2.rds,
##            grf_preds_au_study2.rds,
##            grf_calibration_us_study2.rds,
##            grf_preds_us_study2.rds,
##            Figures/fig_s16.pdf,
##            Figures/fig_s17.pdf,
##            Tables/table_s34.tex,
##            Tables/table_s35.tex
## ============================================================================

source("00_setup.R")

## ----- Data setup -----

mine_dat <-
  combined_dat %>%
  select(
    sample,
    admin_response_id,
    y_trust_m_n,
    y_reputation_m_n,
    y_effective_m_n,
    y_sincerity_m_n,
    Z_vignette_a,
    Z_vignette_order
  ) %>%
  group_by(sample) %>%
  mutate(
    y_mine_aidx = c(
      y_trust_m_n + y_reputation_m_n + y_effective_m_n +
        y_sincerity_m_n
    ),
    y_mine_aidx_s = glass_delta(Y = y_mine_aidx, Z = Z_vignette_a,
                                reference = "Control"),
    Z_vig_type = "Mining"
  ) %>%
  rename(y_aidx = y_mine_aidx,
         y_aidx_s = y_mine_aidx_s,
         Z = Z_vignette_a,
         Z_vig_order = Z_vignette_order,
         y_trust_n = y_trust_m_n,
         y_reputation_n = y_reputation_m_n,
         y_effective_n = y_effective_m_n,
         y_sincerity_n = y_sincerity_m_n)

ent_dat <-
  combined_dat %>%
  select(
    sample,
    admin_response_id,
    y_trust_e_n,
    y_reputation_e_n,
    y_effective_e_n,
    y_sincerity_e_n,
    Z_vignette_b,
    Z_vignette_order
  ) %>%
  group_by(sample) %>%
  mutate(
    y_ent_aidx = c(
      y_trust_e_n + y_reputation_e_n + y_effective_e_n +
        y_sincerity_e_n
    ),
    y_ent_aidx_s = glass_delta(Y = y_ent_aidx, Z = Z_vignette_b,
                               reference = "Control"),
    Z_vig_type = "Entertainment"
  ) %>%
  rename(y_aidx = y_ent_aidx,
         y_aidx_s = y_ent_aidx_s,
         Z = Z_vignette_b,
         Z_vig_order = Z_vignette_order,
         y_trust_n = y_trust_e_n,
         y_reputation_n = y_reputation_e_n,
         y_effective_n = y_effective_e_n,
         y_sincerity_n = y_sincerity_e_n)

stack_dat <-
  bind_rows(mine_dat, ent_dat)

#### ============== Estimate causal forests: Australia ============== ####

au_covs <-
  c(
    "x_edu_college",
    "x_employ_cat",
    "x_gender_f",
    "x_age_cat",
    "x_hhi",
    "x_native",
    "x_pid_comb2",
    "x_state_comb",
    "x_ideo_cat"
  )

set.seed(123)

#### ---------- Substantive v. Control -----

est_sub_dat <-
  stack_dat %>%
  ungroup() %>%
  filter(sample == "AU") %>%
  left_join(combined_dat %>%
              filter(sample == "AU") %>%
              select(admin_response_id, all_of(au_covs)),
            by = "admin_response_id") %>%
  mutate(admin_response_id = factor(admin_response_id)) %>%
  filter(Z != "Symbolic") %>%
  select(admin_response_id,
         all_of(au_covs),
         Z,
         Z_vig_type,
         Z_vig_order,
         y_aidx_s)

index <- complete.cases(est_sub_dat)
X <- est_sub_dat[index, c(au_covs, "Z_vig_type", "Z_vig_order")]

dummies <- model.matrix(~ x_age_cat - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_hhi - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_ideo_cat - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_pid_comb2 - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_employ_cat - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_state_comb - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ Z_vig_type - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ Z_vig_order - 1, data = X)
X <- cbind(X, dummies)

nums <- unlist(lapply(X, is.numeric))
X <- as.matrix(X[, nums])

W <- as.numeric(est_sub_dat$Z[index] == "Substantive")
Y_Obs <- est_sub_dat$y_aidx_s[index]

tau_sub <-
  causal_forest(
    Y = Y_Obs,
    X = X,
    W = W,
    W.hat = 0.5,
    honesty = TRUE,
    num.trees = 5000,
    seed = 123,
    clusters = est_sub_dat$admin_response_id[index]
  )

## Variable importance (used for in-text estimates)
varimp <- c(variable_importance(tau_sub))
names(varimp) <- colnames(X)
ranked_vars <- sort(varimp, decreasing = TRUE)

tau_hat_sub <- predict(tau_sub, estimate.variance = TRUE)
sigma_hat_sub <- sqrt(tau_hat_sub$variance.estimates)

grf_sub <-
  cbind(data.frame(
    tau_hat = tau_hat_sub$predictions,
    ci_upper = tau_hat_sub$predictions + 1.96 * sigma_hat_sub,
    ci_lower = tau_hat_sub$predictions - 1.96 * sigma_hat_sub),
    X
  ) %>%
  mutate(contrast = "Substantive v. Control")

#### ---------- Symbolic v. Control -----

est_sym_dat <-
  stack_dat %>%
  ungroup() %>%
  filter(sample == "AU") %>%
  left_join(combined_dat %>%
              filter(sample == "AU") %>%
              select(admin_response_id, all_of(au_covs)),
            by = "admin_response_id") %>%
  mutate(admin_response_id = factor(admin_response_id)) %>%
  filter(Z != "Substantive") %>%
  select(admin_response_id,
         all_of(au_covs),
         Z,
         Z_vig_type,
         Z_vig_order,
         y_aidx_s)

index <- complete.cases(est_sym_dat)
X <- est_sym_dat[index, c(au_covs, "Z_vig_type", "Z_vig_order")]

dummies <- model.matrix(~ x_age_cat - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_hhi - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_ideo_cat - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_pid_comb2 - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_employ_cat - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_state_comb - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ Z_vig_type - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ Z_vig_order - 1, data = X)
X <- cbind(X, dummies)

nums <- unlist(lapply(X, is.numeric))
X <- as.matrix(X[, nums])

W <- as.numeric(est_sym_dat$Z[index] == "Symbolic")
Y_Obs <- est_sym_dat$y_aidx_s[index]

tau_sym <-
  causal_forest(
    Y = Y_Obs,
    X = X,
    W = W,
    W.hat = 0.5,
    honesty = TRUE,
    num.trees = 5000,
    seed = 123,
    clusters = est_sym_dat$admin_response_id[index]
  )

tau_hat_sym <- predict(tau_sym, estimate.variance = TRUE)
sigma_hat_sym <- sqrt(tau_hat_sym$variance.estimates)

grf_sym <-
  cbind(data.frame(
    tau_hat = tau_hat_sym$predictions,
    ci_upper = tau_hat_sym$predictions + 1.96 * sigma_hat_sym,
    ci_lower = tau_hat_sym$predictions - 1.96 * sigma_hat_sym),
    X
  ) %>%
  mutate(contrast = "Symbolic v. Control")

#### ------ Stack results and export -----

au_forest_test_df <-
  bind_rows(
    "Symbolic v. Control" = broom::tidy(test_calibration(tau_sym)),
    "Substantive v. Control" = broom::tidy(test_calibration(tau_sub)),
    .id = "contrast"
  ) %>%
  mutate(contrast = factor(
    contrast,
    levels = c("Symbolic v. Control", "Substantive v. Control")
  )) %>%
  mutate(across(where(is.numeric), \(x) round(x, 2)))

saveRDS(au_forest_test_df, "grf_calibration_au_study2.rds")

grf_preds_au <-
  bind_rows(grf_sym, grf_sub) %>%
  mutate(contrast = factor(
    contrast,
    levels = c("Symbolic v. Control", "Substantive v. Control")
  )) %>%
  group_by(contrast) %>%
  mutate(order = rank(tau_hat, ties.method = "first")) %>%
  ungroup()

saveRDS(grf_preds_au, "grf_preds_au_study2.rds")

#### ============== Estimate causal forests: United States ============== ####

us_covs <-
  c(
    "x_edu_college",
    "x_employ_cat",
    "x_gender_f",
    "x_age_cat",
    "x_hhi",
    "x_native",
    "x_pid_comb2",
    "x_region",
    "x_ethnicity"
  )

set.seed(123)

#### ---------- Substantive v. Control -----

est_sub_dat <-
  stack_dat %>%
  ungroup() %>%
  filter(sample == "US") %>%
  left_join(combined_dat %>%
              filter(sample == "US") %>%
              select(admin_response_id, all_of(us_covs)),
            by = "admin_response_id") %>%
  mutate(admin_response_id = factor(admin_response_id)) %>%
  filter(Z != "Symbolic") %>%
  select(admin_response_id,
         all_of(us_covs),
         Z,
         Z_vig_type,
         Z_vig_order,
         y_aidx_s)

index <- complete.cases(est_sub_dat)
X <- est_sub_dat[index, c(us_covs, "Z_vig_type", "Z_vig_order")]

dummies <- model.matrix(~ x_age_cat - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_hhi - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_ethnicity - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_employ_cat - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_region - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_pid_comb2 - 1, data = X)
X <- cbind(X, dummies)

nums <- unlist(lapply(X, is.numeric))
X <- as.matrix(X[, nums])

W <- as.numeric(est_sub_dat$Z[index] == "Substantive")
Y_Obs <- est_sub_dat$y_aidx_s[index]

tau_sub <-
  causal_forest(
    Y = Y_Obs,
    X = X,
    W = W,
    W.hat = 0.5,
    honesty = TRUE,
    num.trees = 5000,
    seed = 123,
    clusters = est_sub_dat$admin_response_id[index]
  )

tau_hat_sub <- predict(tau_sub, estimate.variance = TRUE)
sigma_hat_sub <- sqrt(tau_hat_sub$variance.estimates)

grf_sub <-
  cbind(data.frame(
    tau_hat = tau_hat_sub$predictions,
    ci_upper = tau_hat_sub$predictions + 1.96 * sigma_hat_sub,
    ci_lower = tau_hat_sub$predictions - 1.96 * sigma_hat_sub),
    X
  ) %>%
  mutate(contrast = "Substantive v. Control")

#### ---------- Symbolic v. Control -----

est_sym_dat <-
  stack_dat %>%
  ungroup() %>%
  filter(sample == "US") %>%
  left_join(combined_dat %>%
              filter(sample == "US") %>%
              select(admin_response_id, all_of(us_covs)),
            by = "admin_response_id") %>%
  mutate(admin_response_id = factor(admin_response_id)) %>%
  filter(Z != "Substantive") %>%
  select(admin_response_id,
         all_of(us_covs),
         Z,
         Z_vig_type,
         Z_vig_order,
         y_aidx_s)

index <- complete.cases(est_sym_dat)
X <- est_sym_dat[index, c(us_covs, "Z_vig_type", "Z_vig_order")]

dummies <- model.matrix(~ x_age_cat - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_hhi - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_ethnicity - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_employ_cat - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_region - 1, data = X)
X <- cbind(X, dummies)

dummies <- model.matrix(~ x_pid_comb2 - 1, data = X)
X <- cbind(X, dummies)

nums <- unlist(lapply(X, is.numeric))
X <- as.matrix(X[, nums])

W <- as.numeric(est_sym_dat$Z[index] == "Symbolic")
Y_Obs <- est_sym_dat$y_aidx_s[index]

tau_sym <-
  causal_forest(
    Y = Y_Obs,
    X = X,
    W = W,
    W.hat = 0.5,
    honesty = TRUE,
    num.trees = 5000,
    seed = 123,
    clusters = est_sym_dat$admin_response_id[index]
  )

tau_hat_sym <- predict(tau_sym, estimate.variance = TRUE)
sigma_hat_sym <- sqrt(tau_hat_sym$variance.estimates)

grf_sym <-
  cbind(data.frame(
    tau_hat = tau_hat_sym$predictions,
    ci_upper = tau_hat_sym$predictions + 1.96 * sigma_hat_sym,
    ci_lower = tau_hat_sym$predictions - 1.96 * sigma_hat_sym),
    X
  ) %>%
  mutate(contrast = "Symbolic v. Control")

#### ------ Stack results and export -----

us_forest_test_df <-
  bind_rows(
    "Symbolic v. Control" = broom::tidy(test_calibration(tau_sym)),
    "Substantive v. Control" = broom::tidy(test_calibration(tau_sub)),
    .id = "contrast"
  ) %>%
  mutate(contrast = factor(
    contrast,
    levels = c("Symbolic v. Control", "Substantive v. Control")
  )) %>%
  mutate(across(where(is.numeric), \(x) round(x, 2)))

saveRDS(us_forest_test_df, "grf_calibration_us_study2.rds")

grf_preds_us <-
  bind_rows(grf_sym, grf_sub) %>%
  mutate(contrast = factor(
    contrast,
    levels = c("Symbolic v. Control", "Substantive v. Control")
  )) %>%
  group_by(contrast) %>%
  mutate(order = rank(tau_hat, ties.method = "first")) %>%
  ungroup()

saveRDS(grf_preds_us, "grf_preds_us_study2.rds")

#### ============== Generate Fig. S16: Australia ============== ####

grf_preds_au <- read_rds("grf_preds_au_study2.rds")

pred_au <-
  grf_preds_au %>%
  ggplot(., aes(x = tau_hat, y = order)) +
  facet_row(~ contrast) +
  geom_errorbar(aes(xmin = ci_lower, xmax = ci_upper),
                col = "dimgrey",
                alpha = 0.2,
                width = 0,
                linewidth = 0.5,
                orientation = "y") +
  geom_point(size = .5, pch = 20,
             alpha = 0.8,
             fill = "black",
             col = "black") +
  geom_vline(xintercept = 0, linewidth = 0.5, lty = 2) +
  labs(title = "", x = "", y = "") +
  scale_y_continuous(expand = c(0.01, 0.01)) +
  scale_x_continuous(expand = c(0.05, 0.05)) +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(linewidth = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(linewidth = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_blank(),
    axis.line.y = element_blank(),
    panel.border = element_rect(fill = NA),
    strip.background = element_rect(fill = "lightgrey"),
    plot.margin = unit(c(t = -0.5, r = 0.5, b = -0.5, l = -0.5), "lines")
  )

ggsave(
  pred_au,
  filename = paste0(fig_dir, "fig_s16.pdf"),
  width = 5,
  height = 4
)

#### ============== Generate Fig. S17: United States ============== ####

grf_preds_us <- read_rds("grf_preds_us_study2.rds")

pred_us <-
  grf_preds_us %>%
  ggplot(., aes(x = tau_hat, y = order)) +
  facet_row(~ contrast) +
  geom_errorbar(aes(xmin = ci_lower, xmax = ci_upper),
                col = "dimgrey",
                alpha = 0.2,
                width = 0,
                linewidth = 0.5,
                orientation = "y") +
  geom_point(size = .5, pch = 20,
             alpha = 0.8,
             fill = "black",
             col = "black") +
  geom_vline(xintercept = 0, linewidth = 0.5, lty = 2) +
  labs(title = "", x = "", y = "") +
  scale_y_continuous(expand = c(0.01, 0.01)) +
  scale_x_continuous(expand = c(0.05, 0.05)) +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(linewidth = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(linewidth = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_blank(),
    axis.line.y = element_blank(),
    panel.border = element_rect(fill = NA),
    strip.background = element_rect(fill = "lightgrey"),
    plot.margin = unit(c(t = -0.5, r = 0.5, b = -0.5, l = -0.5), "lines")
  )

ggsave(
  pred_us,
  filename = paste0(fig_dir, "fig_s17.pdf"),
  width = 5,
  height = 4
)

#### ============== Table S34: Omnibus tests for heterogeneity ============== ####

grf_calibrate_au <- read_rds("grf_calibration_au_study2.rds")
grf_calibrate_us <- read_rds("grf_calibration_us_study2.rds")

au_rows <- grf_calibrate_au %>%
  filter(term == "differential.forest.prediction") %>%
  arrange(contrast)

us_rows <- grf_calibrate_us %>%
  filter(term == "differential.forest.prediction") %>%
  arrange(contrast)

tex_row <- function(label, est, se, stat, pval) {
  paste0("\\quad ", label, " & ", est, " & ", se, " & ", stat, " & ", pval, " \\\\")
}

writeLines(
  c(
    "\\multicolumn{5}{l}{\\textit{Australia:}} \\\\",
    tex_row(as.character(au_rows$contrast), au_rows$estimate,
            au_rows$std.error, au_rows$statistic, au_rows$p.value),
    "\\multicolumn{5}{l}{\\textit{United States:}} \\\\",
    tex_row(as.character(us_rows$contrast), us_rows$estimate,
            us_rows$std.error, us_rows$statistic, us_rows$p.value)
  ),
  paste0(tab_dir, "table_s34.tex")
)

#### ============== In-text estimates: S2.3.1 ============== ####

## Variable importance percentages: AU Substantive v. Control
varimp_pct <- round(100 * varimp, 0)
cat("Variable importance (% of splits) for top 5 covariates:\n")
print(sort(varimp_pct[names(ranked_vars[1:5])], decreasing = TRUE))

## Subgroup proportions: unique AU respondents
au_unique <- combined_dat %>%
  filter(sample == "AU") %>%
  distinct(admin_response_id, .keep_all = TRUE)

cat("% Left partisans:",
    round(100 * mean(au_unique$x_pid_comb2 == "Left", na.rm = TRUE), 0), "\n")
cat("% Left ideology:",
    round(100 * mean(au_unique$x_ideo_cat == "Left", na.rm = TRUE), 0), "\n")
cat("% Moderate ideology:",
    round(100 * mean(au_unique$x_ideo_cat == "Moderate", na.rm = TRUE), 0), "\n")
cat("% Right ideology:",
    round(100 * mean(au_unique$x_ideo_cat == "Right", na.rm = TRUE), 0), "\n")
cat("% College graduates:",
    round(100 * mean(au_unique$x_edu_college == 1, na.rm = TRUE), 0), "\n")

#### ============== Table S35: AU CATE subgroups ============== ####

## Reconstruct AU Substantive dataset with all covariates for subgroup 
## regressions
est_sub_dat <-
  stack_dat %>%
  ungroup() %>%
  filter(sample == "AU") %>%
  left_join(combined_dat %>%
              filter(sample == "AU"),
            by = "admin_response_id") %>%
  mutate(admin_response_id = factor(admin_response_id)) %>%
  filter(Z != "Symbolic")

cate_tab <-
  lm_robust(
    y_aidx_s ~ Z + Z_vig_order + Z_vig_type,
    data = est_sub_dat %>% filter(x_pid_comb2 == "Left"),
    clusters = admin_response_id
  ) %>%
  tidy() %>%
  filter(term == "ZSubstantive") %>%
  mutate(across(c(estimate, std.error), \(x) round(x, 2))) %>%
  mutate(X = "Partisanship: Left", estimand = "CATE1") %>%
  bind_rows(
    lm_robust(
      y_aidx_s ~ Z + Z_vig_order + Z_vig_type,
      data = est_sub_dat %>% filter(x_pid_comb2 != "Left"),
      clusters = admin_response_id
    ) %>%
      tidy() %>%
      filter(term == "ZSubstantive") %>%
      mutate(across(c(estimate, std.error), \(x) round(x, 2))) %>%
      mutate(X = "Partisanship: Independent/Right", estimand = "CATE2")
  ) %>%
  bind_rows(
    lm_robust(
      y_aidx_s ~ Z * X + Z_vig_order + Z_vig_type,
      data = est_sub_dat %>% mutate(X = as.numeric(x_pid_comb2 == "Left")),
      clusters = admin_response_id
    ) %>%
      tidy() %>%
      filter(term == "ZSubstantive:X") %>%
      mutate(across(c(estimate, std.error), \(x) round(x, 2))) %>%
      mutate(X = "Partisanship: Left", estimand = "Interaction")
  ) %>%
  bind_rows(
    lm_robust(
      y_aidx_s ~ Z + Z_vig_order + Z_vig_type,
      data = est_sub_dat %>% filter(x_ideo_cat == "Left"),
      clusters = admin_response_id
    ) %>%
      tidy() %>%
      filter(term == "ZSubstantive") %>%
      mutate(across(c(estimate, std.error), \(x) round(x, 2))) %>%
      mutate(X = "Ideology: Left", estimand = "CATE1")
  ) %>%
  bind_rows(
    lm_robust(
      y_aidx_s ~ Z + Z_vig_order + Z_vig_type,
      data = est_sub_dat %>% filter(x_ideo_cat %in% c("Moderate", "Right")),
      clusters = admin_response_id
    ) %>%
      tidy() %>%
      filter(term == "ZSubstantive") %>%
      mutate(across(c(estimate, std.error), \(x) round(x, 2))) %>%
      mutate(X = "Ideology: Moderate/Right", estimand = "CATE2")
  ) %>%
  bind_rows(
    lm_robust(
      y_aidx_s ~ Z * X + Z_vig_order + Z_vig_type,
      data = est_sub_dat %>% mutate(X = as.numeric(x_ideo_cat == "Left")),
      clusters = admin_response_id
    ) %>%
      tidy() %>%
      filter(term == "ZSubstantive:X") %>%
      mutate(across(c(estimate, std.error), \(x) round(x, 2))) %>%
      mutate(X = "Ideology: Left", estimand = "Interaction")
  ) %>%
  bind_rows(
    lm_robust(
      y_aidx_s ~ Z + Z_vig_order + Z_vig_type,
      data = est_sub_dat %>% filter(x_edu_college == 1),
      clusters = admin_response_id
    ) %>%
      tidy() %>%
      filter(term == "ZSubstantive") %>%
      mutate(across(c(estimate, std.error), \(x) round(x, 2))) %>%
      mutate(X = "College: Yes", estimand = "CATE1")
  ) %>%
  bind_rows(
    lm_robust(
      y_aidx_s ~ Z + Z_vig_order + Z_vig_type,
      data = est_sub_dat %>% filter(x_edu_college == 0),
      clusters = admin_response_id
    ) %>%
      tidy() %>%
      filter(term == "ZSubstantive") %>%
      mutate(across(c(estimate, std.error), \(x) round(x, 2))) %>%
      mutate(X = "College: No", estimand = "CATE2")
  ) %>%
  bind_rows(
    lm_robust(
      y_aidx_s ~ Z * X + Z_vig_order + Z_vig_type,
      data = est_sub_dat %>% mutate(X = as.numeric(x_edu_college)),
      clusters = admin_response_id
    ) %>%
      tidy() %>%
      filter(term == "ZSubstantive:X") %>%
      mutate(across(c(estimate, std.error), \(x) round(x, 2))) %>%
      mutate(X = "College: Yes", estimand = "Interaction")
  ) %>%
  bind_rows(
    lm_robust(
      y_aidx_s ~ Z + Z_vig_order + Z_vig_type,
      data = est_sub_dat %>% filter(x_ideo_cat == "Moderate"),
      clusters = admin_response_id
    ) %>%
      tidy() %>%
      filter(term == "ZSubstantive") %>%
      mutate(across(c(estimate, std.error), \(x) round(x, 2))) %>%
      mutate(X = "Ideology: Moderate", estimand = "CATE1")
  ) %>%
  bind_rows(
    lm_robust(
      y_aidx_s ~ Z + Z_vig_order + Z_vig_type,
      data = est_sub_dat %>% filter(x_ideo_cat != "Moderate"),
      clusters = admin_response_id
    ) %>%
      tidy() %>%
      filter(term == "ZSubstantive") %>%
      mutate(across(c(estimate, std.error), \(x) round(x, 2))) %>%
      mutate(X = "Ideology: Left/Right", estimand = "CATE2")
  ) %>%
  bind_rows(
    lm_robust(
      y_aidx_s ~ Z * X + Z_vig_order + Z_vig_type,
      data = est_sub_dat %>% mutate(X = as.numeric(x_ideo_cat == "Moderate")),
      clusters = admin_response_id
    ) %>%
      tidy() %>%
      filter(term == "ZSubstantive:X") %>%
      mutate(across(c(estimate, std.error), \(x) round(x, 2))) %>%
      mutate(X = "Ideology: Moderate", estimand = "Interaction")
  ) %>%
  bind_rows(
    lm_robust(
      y_aidx_s ~ Z + Z_vig_order + Z_vig_type,
      data = est_sub_dat %>% filter(x_pid_comb2 == "Right"),
      clusters = admin_response_id
    ) %>%
      tidy() %>%
      filter(term == "ZSubstantive") %>%
      mutate(across(c(estimate, std.error), \(x) round(x, 2))) %>%
      mutate(X = "Partisanship: Right", estimand = "CATE1")
  ) %>%
  bind_rows(
    lm_robust(
      y_aidx_s ~ Z + Z_vig_order + Z_vig_type,
      data = est_sub_dat %>% filter(x_pid_comb2 != "Right"),
      clusters = admin_response_id
    ) %>%
      tidy() %>%
      filter(term == "ZSubstantive") %>%
      mutate(across(c(estimate, std.error), \(x) round(x, 2))) %>%
      mutate(X = "Partisanship: Left/Independent", estimand = "CATE2")
  ) %>%
  bind_rows(
    lm_robust(
      y_aidx_s ~ Z * X + Z_vig_order + Z_vig_type,
      data = est_sub_dat %>% mutate(X = as.numeric(x_pid_comb2 == "Right")),
      clusters = admin_response_id
    ) %>%
      tidy() %>%
      filter(term == "ZSubstantive:X") %>%
      mutate(across(c(estimate, std.error), \(x) round(x, 2))) %>%
      mutate(X = "Partisanship: Right", estimand = "Interaction")
  ) %>%
  select(estimand, X, estimate, std.error, statistic, p.value)

cate_tab %>%
  mutate(
    group = c(
      rep("Left partisanship:", 3),
      rep("Left ideology:", 3),
      rep("College graduate:", 3),
      rep("Moderate ideology:", 3),
      rep("Right Partisanship:", 3)
    ),
    entry = table_entry(est = estimate, se = std.error),
    pval = case_when(
      p.value < 0.001 ~ "< 0.001",
      round(p.value, 3) == 0.001 ~ "0.001",
      p.value > 0.99 ~ "> 0.99",
      .default = sprintf("%.3f", p.value)
    )
  ) %>%
  select(estimand, group, entry, pval) %>%
  pivot_wider(names_from = estimand, values_from = c(entry, pval)) %>%
  select(
    group,
    entry_CATE1, pval_CATE1,
    entry_CATE2, pval_CATE2,
    entry_Interaction, pval_Interaction
  ) %>%
  xtable(.) %>%
  print(
    include.rownames = FALSE,
    include.colnames = FALSE,
    hline.after = c(),
    only.contents = TRUE,
    type = "latex",
    sanitize.text.function = identity,
    comment = FALSE,
    file = paste0(tab_dir, "table_s35.tex")
  )

## ---- Session Info -----------------------------------------------------------
sessionInfo()
